Hybrid density functional calculations of the band gap of Ga x Ini_ x N 
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Recent theoretical work has provided evidence that hybrid functionals, which include a fraction 
of exact (Hartree Fock) exchange in the density functional theory (DFT) exchange and correlation 
terms, significantly improve the description of band gaps of semiconductors compared with local and 
semilocal approximations. Based on a recently developed order- N method for calculating the exact 
exchange in extended insulating systems, we have implemented an efficient scheme to determine 
the hybrid functional band gap. We use this scheme to study the band gap and other electronic 
properties of the ternary compound Ini_ x Ga x N using a 64-atom supercell model. 
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The design of novel functional semiconductors with 
given values of the energy band gap is an area of in- 
tense researchi* 2 -^^. In particular, much attention is 
focused on the band gap engineering of group-Ill nitride 
semiconductors, whose remarkable optical properties are 
important for optoelectronic device applications 6,7 . To 
guide the search for compounds with tailored properties 1 , 
experimental studies are often accompanied by electronic 
structure calculations based on density functional theory 
(DFT) 8 . For these calculations, the local-density (LDA) 
or generalized gradient approximations(GGA) are typi- 
cally used. Due to the derealization error of the LDA 
and GGA exchange and correlation functionals, however, 
these approaches severely underestimate the materials 
band gapsSaiS. 

As shown by several recent studies 11 , a significant im- 
provement in the description of semiconductor and insu- 
lator band gaps is generally obtained by using hybrid 
functionals 12 , in which some exact (Hartree- Fock) ex- 
change is mixed into the exchange and correlation func- 
tional. This reduces the derealization and derivative dis- 
continuity errors of (semi) local functionals 2 * ^ 1Q i 11 . How- 
ever, because of the considerable computational cost of 
evaluating the non-local exact exchange term, hybrid 
functionals have been mostly applied to systems with 
small unit cells 13 . For the modeling of systems where a 
large supercell is needed, an additional screened exchange 
approximation is usually made to relieve the computa- 
tional burden 2 -^. 

Recently Wu et al (WSC)^ introduced an order-TV 
method to calculate the exact exchange in extended in- 
sulating systems. The WSC method is based on a lo- 
calized Wannier function representation of the occupied 
(valence) space, so that the exchange interaction between 
two orbitals decays rapidly with the distance between 
their centers. A truncation can thus be introduced, which 
greatly reduces the computational cost. The effectiveness 
of the WSC method was demonstrated by ground state 
electronic minimizations for crystalline silicon in super- 
cells with 64 and 216 atoms. 



In this paper, we extend the WSC scheme to compute 
hybrid functional band gaps. To this end, the system's 
first (few) empty conduction state (s) is (are) determined 
starting from the ground state calculated via the WSC 
method. With hybrid functionals, this requires the com- 
putation of the pair exchange between the empty state 
and each valence orbital. Even though the empty state 
is delocalized, the product between this state and a va- 
lence orbital is well localized, so that the corresponding 
exchange interaction can be truncated as in the origi- 
nal WSC methodM. We apply our scheme to determine 
the band gap of Ini-^Ga^N, a ternary nitride semicon- 
ductor of great technological interest, and of its par- 
ent compounds, InN and GaN, using the PBEO hybrid 
functional^. Our results show that, compared to the 
semi-local PBE functional, PBEO gives a considerably 
improved description of the band gap, as well as of the 
cation d state binding energy, which is also poorly de- 
cribed by the semilocal functionals. 

The PBEO hybrid functional is constructed by mixing 
25% of exact exchange with the GGA-PBE exchange 16 , 
while the correlation potential is still represented by the 
corresponding functional in PBEii, 
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Here E x denotes the exact exchange energy, E PBE is 
the PBE exchange, and E PBE is the PBE correlation 
functional. E x has the usual Hartree-Fock form in terms 
of one-electron orbitals. In the WSC method, this term 
is expressed in terms of localized Wannier orbitals {(fi}. 
These are obtained through an unitary transformation 
of the delocalized Bloch states {ipi} corresponding to oc- 
cupied bands. In particular, we use maximally localized 
Wannier functions (MLWFs)^, which are exponentially 
localized. In this way, a significant truncation in both 
number and size of exchange pairs can be achieved in 
real space. 

We now turn to the calculation of the band gap. In 
extended insulating systems the band gap is simply given 
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by the difference between the eigenvalue of the highest 
occupied and the lowest empty state. Once the ground 
state has been minimized self-consistently, the eigenvalue 
of the empty state cp e can be obtained through a simple 
non-selfconsistent calculation. With the hybrid PBEO 
functional, the equation for (p e is 

( - ^V 2 + F ion (r) +^ H [p val (r)] + ^ PBE [p val (r)] 
+ K PBE [p val W]) x^(r) + i J Vr\ry)Mr f )dr f 

= £ e ^e0), (2) 

In the above expression we have assumed, for simplic- 
ity, a closed-shell system with N/2 doubly occupied one- 
electron states (extension to spin-polarized systems is 
straightforward); Vr and VJ on are the Hartree and the 
ionic (pseudo-)potentials, respectively; V BBE and V BBE 
are the PBE exchange and correlation potentials. We 
note that Vn, ^ PBE an d V BBE are fixed operators as 
they only depend on the (fixed) valence charge density 
p val (r) = Y^ CC <Pj ( r )^i( r )- Finally, the non-local exact 
exchange potential 1/J al (r, r') is given by: 

^(r,0 = -2E i r ! i > (3) 
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where the sum runs over all the occupied states. This 
potential describes the exchange interaction between the 
empty state and each of the valence MLWFs {fij}. 

The action of ^ al (r, r') on the empty state cp e in 
Eq. ([3]) is given by: 

3 
occ 

= -2j> ei (r)&(r) (4) 
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Here v e j is the Coulomb potential originating from the 
"exchange charge" p e j = lpj(r')(p e (r'), and satisfies the 
Poisson equation: 

V 2 v ej = -47rp e j (5) 

It is important to note that, while the empty eigenstate 
of Eq. ([2]) is Bloch like and delocalized in real space, the 
exchange pair density p e j is confined by the valence ML- 
WFs that are well localized in real space. As a result, the 
Poisson equation, Eq. ([5]), and the action of the exchange 
operator, Eq. (jH), need only be solved in the region where 

We have implemented the above computational pro- 
cedure for calculating the PBEO band gap in the CP 
code of the Quantum-ESPRESSO packaged The pro- 
cedure works as a post processing feature following a 
PBEO ground state calculation by the MLWF-based 
WSC method. In this work, we use it to calculate the 
electronic structure, particularly the band gap, of GaN, 



InN, and Ini-^Ga^N in the zincblende phase. These sys- 
tems are computationally challenging because InN and 
In-rich Ini-^Ga^N are incorrectly predicted to be metal- 
lic by standard GGA calculations. 

The calculations were performed using a 64-atom cubic 
supercell to model both Ini-^Ga^N and its parent com- 
pounds, GaN and InN. For each Ga concentration x in 
the ternary Ini-^Ga^N compound, only a few selected 
atomic configurations were considered, with no specific 
treatment of disorder effects, as e.g. in Refs. 00; within 
our limited sampling, a very weak dependence of the cal- 
culated band gap on the specific cation arrangement was 
observed. For direct comparison with experiments and 
other theoretical results, the experimental lattice con- 
stants of GaN (a = 4.50 A) and InN (a = 4.98 A) were 
used, while the lattice parameter of the alloy was deter- 
mined by linear interpolation. 

TABLE I: Pseudopotential generation parameters. Here 
"ref." refers to the reference state occupation, r c refers to 
the cut-off radius, q c is the cut-off wavevector and Nb is the 
number of Bessel functions used for each channel (see Ref. l28h . 
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ref. 


2.0 


3.0 






r c 


1.30 


1.30 






Qc 


7.50 


7.50 






N B 


10 


10 




Ga 


ref. 


2.0 


1.0 


10.0 




r c 


1.80 


2.20 


1.80 




Qc 


8.00 


8.00 


8.36 
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10 


In 


ref. 


2.0 


1.0 


10.0 




r c 


1.90 


2.30 


1.80 




Qc 


8.00 


8.00 


8.00 




N B 


8 


8 


8 



Table [J shows the reference states and cut-off radii 
used to construct the pseudopotentials used in this study. 
All pseudopotentials were generated using the OPIUM 
code^i. 

Unlike with traditional density functional theory, 
Hartree- Fock pseudopotentials require extra care in their 
construction. This arises from the non-local form of the 
Hartree- Fock exchange potential 22 ,23 ! 24 ! 25 . The presence 
of the non-local exchange potential in Hartree-Fock or 
Hartree- Fock/DFT hybrids will often yield pseudopoten- 
tials with an unphysical, long-range tail. A correction 
procedure is necessary to remove this tail and restore 
the correct long-range behavior of the pseudopotential 
while maintaining the eigenvalue spectrum and logarith- 
mic derivatives. Recent work 2 2 i 26 i 27 has shown that this 
approach yields highly accurate Hartre-Fock pseudopo- 
tentials. 

The pseudopotentials were norm-conserving/RRKJ 
type 28 and were generated from self-consistent PBEO all- 
electron reference states using the approach of Ref. [2?]. 
The Ga and In pseudopotentials were obtained from 
scalar-relativistic solutions, while the N pseudopotential 
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was non-relativistic. The local potential was the s chan- 
nel for all cases. The semi-core d electrons were treated 
as valence electrons in In and Ga (this corresponds to 576 
valence electrons, i.e. 288 occupied states, in the 64-atom 
supercell). The plane- wave energy cutoff was 70 Ry and 
the Brillouin zone was sampled at the Y point. Atomic 
positions in the supercell were relaxed at the GGA-PBE 
level. 

The PBEO ground state was determined by the WSC 
method, using MLWFs to calculate the exchange inter- 
action among valence electrons^. While the MLWFs 
generated from the PBE ground state often give an ex- 
cellent initial guess for the PBEO calculations, for InN 
and In rich Ga^Ini^N alloy configurations, the PBE 
ground state shows an incorrect ordering of the energy 
bands. For this reason, instead of PBE Wannier orbitals 
we used a set of fictitious localized orbitals at the guess 
bonding centers as the trial solutions for Eq. ([2]). This 
procedure was essential to obtain the PBEO ground state 
with correct symmetry for InN and In rich Ga^Ini-^N. 
In the empty state calculations, for each PBEO ground 
state MLWF we first defined an orthorhombic box such 
that outside this box p e j (r) is smaller than a given cut-off 
value p cut ; we take this cut-off equal to 2 x 10 -4 bohr -3 in 
the present work. Then Eq. ([5]) is solved by the conjugate 
gradient method 1 -, and for each pair p e j formed by the 
empty state and a PBEO ground state MLWF its action 
Eq. ([I]) is applied only inside the above truncated box. 
Finally with this 2}®(r), Eq. ([2]) is solved via a damped 
second order Car-Parrinello dynamics 19 . 




FIG. 1: (Color online.) Isosurfaces of typical d-like and sp 3 
like Wannier orbitals in the InN (on the left) and GaN (on 
the right) 64-atom supercell. The Ga, In and N atoms are 
denoted by the green, red and blue spheres respectively. 

Representative MLWFs for InN in its PBEO ground 
state are shown in Fig. [TJ Two types of valence MLWFs 
are present in our calculations, d-\ike Wannier orbitals 
centered at the In sites, and covalent sp 3 -\ike orbitals 
centered between the cations and the anions. As one 
can see from the figure, the d-\ike orbitals originating 
from the cation semi-core states are more localized than 
the sp 3 -like ones. The valence MLWFs are qualitatively 
similar for GaN, except for a slightly more pronounced 
localization related to the larger band gap. 

The band structure properties of GaN and InN that 
result from our PBEO-MLWFs calculations are summa- 



TABLE II: Valence band width, band gap and average d-band 
binding energy (eV) of GaN and InN. 







VBW 


Eg 


E d 


GaN 


PBEO-MLWFs 


17.70 


3.52 


-16.16 




PBE 


16.14 


1.60 


-13.62 




PBEO, plane waves a 


17.72 


3.61 






GW b 




3.53 


-16.5 




Experiment c 




3.3 


-17.7 


InN 


PBEO-MLWFs 


17.04 


1.09 


-15.30 




PBE results 


15.04 


-0.04 


-13.48 




GW 6 




0.78 


-15.3 




Experiment 




0.61 6 


-16.0 C 



a Reciprocal space method in PWSCF (Ref. 1 2 (J) in 2-atom cell and 
4x4x4 k points 
b Reference [29l 
c Reference Om 



rized in Table HI Here we report the valence band width 
(VBW), the band gap E g and the average d-band bind- 
ing energy Ed, and compare them to PBE calculations 
(performed with the same 64-atom supercell used for the 
PBEO calculations) and experimental results. For further 
comparison, we also report the results of PBEO calcula- 
tions performed using the reciprocal space implementa- 
tion in Ref. [20|; we can see that the agreement between 
these results and our MLWF-based calculations is very 
good. From Table [Til it appears that the GGA-PBE re- 
sults significantly overestimate the energetic position of 
the cation d-bands. Because of the pd repulsion, the 
overestimated d bands level in turn pushes the p band 
upwards, resulting in an underestimated band gap. For 
InN, this effect leads to a wrong ordering of the Ti c and 
Ti5 V energy levels, and thus to the incorrect prediction 
of a metallic ground state. In the PBEO calculations, the 
inclusion of exact exchange reduces the derealization er- 
ror. As shown by Table HI1 the PBEO VBW is larger and 
the (i-bands level shifts downwards, in better agreement 
with the experiment. In turn, this leads to a considerable 
improvement of the band gaps of both InN and GaN with 
respect to experiment; in particular, the PBEO band gap 
becomes 1.09 eV for InN. It is also worth noticing that 
calculation of the PBEO band gap using a PBE pseudopo- 
tential yields a ~ 0.2 eV smaller value than that obtained 
with the PBEO pseudopotential. 

Besides confirming the good performance of hybrid 
functionals for band gap predictions, the above results 
for InN and GaN provide evidence of the reliability of 
our procedure for calculating the PBEO band gap. We 
have thus applied this procedure to the study of the 
ternary Ini-^Ga^N compound, a system for which the 
standard reciprocal space approach to calculate the exact 
exchange would be extremely cumbersome. Instead, our 
order-TV scheme is well suited to treat systems for which 
large supercells are needed. Using a 64-atom supercell, 
we then considered Ini-^Ga^N models with 1(31), 2(30), 
3(29), 4(28), 16(16), 28(4), 29(3), 30(2) and 31(1) Ga(In) 
cations, which correspond to x = 0.031, 0.063, 0.094, 



4 



3.5 
f 3.0 

"D 2.0 
03 

-Q 1.5 
o 

LJJ 

CD 10 
CL 

0.5 

0.0 

6 



^ 5 
> 

CD 4 
CD 

?» 

c 
05 

n 2 
o 

LU 

DQ h 
CL 1 



. PBEO 




„ pop 
■ iDl 




- — Exper. 






* 

w 




■ 

■ 


■ 

T " " " 1 i i i 


(a) 


- VBM 


• 

• 


" . CBM 


• 

• 


• 




■ 







(b) 


0.0 0.2 0.4 0.6 


0.8 1.0 



Ga fraction x in ln-| -xGa x N 

FIG. 2: (Color online.) (a) PBEO, PBE and experimental 
band gap of dependence Ga fraction x (b) Valence band max- 
imum (VBM) and conduction band minimum as a function 
of Ga fraction x in Ini-^GaxN 



0.125, 0.5, 0.875, 0.906, 0.938, and 0.969. For each 
value of x and a given configuration of Ga(In) atoms, 
the atomic positions were relaxed at the PBE level. The 
computed PBEO band gap of Ini-^Ga^N as a function 
of the Ga fraction x is shown in Fig. Efa), together with 
experimental 31 and PBE results. We can see that PBE 
not only significantly underestimates the band gap but 
incorrectly shows a metallic ground state for x < 0.5. By 
contrast, a direct band gap at the T point is found for all 
values of x at the PBEO level. Moreover, PBEO predicts 
a large band gap bowing effect, in qualitative agreement 
with the experiment* 3 ^. The band gap can be fitted to 
the quadratic form 



^alloy 



xE\ 



GaN 



+ (1 - x)E lnN - x(l - x)b (6) 



from which a bowing coefficient fr PBE0 = 1.63 eV can 
be extracted, similar to the value, 1.67 eV, found in 
previous screened-exchange density functional (sx-LDA) 
calculations 4 . However, this is somewhat larger than the 
experimental value fr expt = 1.43 eV— , likely because of 
the overestimated PBEO band gap for the In-rich com- 
pounds. To gain more insight into the origin of the large 
band gap bowing, we have examined how the valence 
band maximum (VBM) and conduction band minimum 
(CBM) depend separately on x, see Fig. [3]^a). In this 
analysis, the average electrostatic potential was taken as 



the reference for the band alignment. It can be seen 
that the VBM increases almost linearly with x, whereas 
the CBM shows a stronger nonlinear increase which is 
responsible for the large bowing coefficient of the alloy. 




FIG. 3: (Color online.) Isosurfaces of PBEO eigenstate (a) 
In3Ga2gN where 3 In atoms forms a zigzag chain structure; 
(b) Ga3ln2gN where 3 Ga atoms forms a zigzag chain N atoms 
are denoted by red, orange and blue spheres respectively 



The electronic states in proximity of the VBM 
are important for the pholuminescence properties of 
Ini-^Ga^N. These states have the character of p orbitals 
localized at the N sites. Previous theoretical studies 
of Ini-^Ga^N found that in Ga-rich alloys the ampli- 
tude of these states is enhanced at N sites close to In 
impurities 4 ^, suggesting a localization of photoexcited 
holes at such sites. This interesting result is confirmed by 
our PBEO hybrid calculations. The enhancement, or hole 
localization, is particularly evident when the In impuri- 
ties are clustered to form a zigzag In-N-In-N-In chain, as 
shown in Fig. [3](a). This localization has been suggested 
to be the reason of the high efficiency of Ini-^Ga^N 
based emitting devices 4,6 . Interestingly, we found that 
there is an opposite effect for the case of Ga impurities 
in In rich alloys. Here, a reduction of the p states at the 
N sites along the Ga-N-Ga-N-Ga-N chain is observed, see 
Fig. lib). 



In conclusion, we have described an efficient procedure 
to calculate the band gap of extended insulating systems 
using hybrid functionals. This procedure is based on 
the recently developed WSC order-TV method, in which 
the Hartree Fock exchange is calculated using MLWFs, 
and can therefore be used to study the band gap and 
other electronic properties of systems with large unit 
cells. We have demonstrated the effectiveness of our ap- 
proach by a study of the band gap of a ternary com- 
pound, Ini-^Ga^N, that we have modeled using a 64- 
atom supercell. Hybrid functional results for this impor- 
tant material are here reported for the first time, with- 
out the approximation of screened exchange, and show 
a much better agreement with experiment than conven- 
tional DFT-GGA or LDA calculations. Our approach 
can be widely used for the band gap engineering problem 
in semiconductor alloys. 
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